An Introduction to Image Segmentation

Image segmentation is the process of partitioning an image into several coherent sub-regions according to some extracted features (e.g., color, brightness, texture attributes). It plays an important role in a broad range of applications, including video surveillance, augmented reality, recognition tasks (face, fingerprints, iris, etc.), autonomous vehicles (navigable surface, object detection, etc.) or medical image analyse, just to name a few. Segmentation aims to extract meaningful information for easier and enhanced analysis; as an illustrative example, a diagnosis based on a set of features extracted from a segmented CT scan of a lung will most probably provide more accurate results than another one where the lung is not isolated.

Many distinct segmentation methods and algorithms have been discussed in the literature, and there are several ways to classify them. For instance, they can be classified into those: semantic segmentation, which performs pixel-level labeling with a set of object categories (e.g., car, human) for all image pixels, and instance segmentation, which extends semantic segmentation scope further by detecting and delineating each object of interest in the image. Also, there is another type called panoptic segmentation which is the unified version of the two basic segmentation processes.

An example of diferent types of image segmentation.

An example of diferent types of image segmentation.

They can also be split into three categories: manual, semiautomatic, and automatic. While an advantage of manual segmentation method is that we can use expert knowledge (as it is usually performed by an expert), its drawbacks are that it is very time consuming and prone to intra and interobserver variability, which can result in a large difference in the extracted features.

Semiautomatic segmentation tries to solve some of the problems related to manual segmentation; by using algorithms, the effort and time spent by the user can be reduced. Semiautomatic algorithms aim to reduce the inter and intraobserver variability. However, interobserver variability will still be present, as the manual part of the segmentation and the settings of the algorithm influence the result.

Automatic segmentation methods do not rely on user interaction, and can be split into two categories: learning and nonlearning-based. Deep learning recently has become very popular as a learning-based method, where the segmentation is performed by a neural network that is trained with labeled examples. The advantage of automatic segmentation is that once the method has been constructed, the segmentations can be performed relatively quickly. Moreover, they produce segmentations that are consistent and reproducible. The disadvantage of deep learning is that it usually requires a lot of labeled data to train an accurate model, as well as a long training time and specialized hardware to construct the model.

According to its architecture, image segmentation algorithms can be grouped into different types, such as:

DICOM Processing and Segmentation in Python

DICOM (Digital Imaging and Communications in Medicine) is a standard protocol for the management and transmission of medical images and related data and is used in many healthcare facilities. This section is dedicated to visualizing Achilles tendon CT images (DICOM files) in 1D (by histogram), 2D, and 3D. Finally, we will create segmentation masks that remove all voxels except for the Achilles tendon. In order to achieve this, we’ll run a Python code based on the one from [2] in RStudio using the library reticulate.

Getting Ready

The reticulate package provides a comprehensive set of tools for interoperability between Python and R. Thus, we’ll start by installing and loading it (see [3] and [4]). Now we are almost ready to run Python code in RStudio. Before doing this, we must install some Python modules (if they have not been installed yet) by running the following command lines (in Terminal):

pip install pydicom

conda install plotly

conda install scikit-image (make sure you are using scikit-image 0.13 or above)

pip install sklearn

pip install matplotlib

Import packages

Once we are working in a Python environment, let us begin by importing the following packages:

import numpy as np
import pydicom
import os
import matplotlib.pyplot as plt
from glob import glob
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import scipy.ndimage
from skimage import morphology
from skimage import measure
from skimage.transform import resize
from sklearn.cluster import KMeans
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from plotly.tools import FigureFactory as FF
from plotly.graph_objs import *
from pydicom.pixel_data_handlers.util import apply_modality_lut
import plotly.figure_factory
import matplotlib

To specify the input and output data folders path:

data_path = "./Images/DICOM/AX DP 1RA68062501 (id=1)/FILESET"
output_path = working_path = "./Images/DICOM/AX DP 1RA68062501 (id=1)"
g = glob(data_path + '/*')
# Print out the first 5 file names to verify we are in the right folder
print ("Total of %d DICOM images.\nFirst 5 filenames:" % len(g))
print ('\n'.join(g[:5]))

Helper Functions

  • load_scan: loads all DICOM images from a folder into a list for manipulation.
  • get_pixels_hu: converts raw voxel and pixel values into Houndsfield Units, which are a relative quantitative measurement of radiodensity used by radiologists in the interpretation of computed tomography (CT) images.

Dispaying an Image Stack

We’ll be skipping every 2 slices to get a representative look at the study.

id = 1
imgs_to_process = np.load(output_path+'/fullimages_{}.npy'.format(id))
def sample_stack(stack, rows=4, cols=5, start_with=0, show_every=2):
    fig,ax = matplotlib.pyplot.subplots(rows,cols,figsize=[12,12])
    #matplotlib.pyplot.subplots_adjust(hspace=0.25)
    for i in range(rows*cols):
        ind = start_with + i*show_every
        ax[int(i/cols),int(i % cols)].set_title('slice %d' % ind)
        ax[int(i/cols),int(i % cols)].imshow(stack[ind],cmap='gray')
        ax[int(i/cols),int(i % cols)].axis('off')
    plt.show()
sample_stack(imgs_to_process)

Resampling

print ("Slice Thickness: %f" % patient[0].SliceThickness)
## Slice Thickness: 2.999637
print ("Pixel Spacing (row, col): (%f, %f) " % (patient[0].PixelSpacing[0], patient[0].PixelSpacing[1]))
## Pixel Spacing (row, col): (0.198900, 0.198900)
id = 1
imgs_to_process = np.load(output_path+'/fullimages_{}.npy'.format(id))
def resample(image, scan, new_spacing=[1,1,1]):
    # Determine current pixel spacing
    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))
    spacing = np.array(list(spacing))
    resize_factor = spacing / new_spacing
    new_real_shape = image.shape * resize_factor
    new_shape = np.round(new_real_shape)
    real_resize_factor = new_shape / image.shape
    new_spacing = spacing / real_resize_factor
    
    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)
    
    return image, new_spacing
print ("Shape before resampling\t", imgs_to_process.shape)
## Shape before resampling   (40, 704, 704)
imgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])
## <string>:11: DeprecationWarning:
## 
## Please use `zoom` from the `scipy.ndimage` namespace, the `scipy.ndimage.interpolation` namespace is deprecated.
print ("Shape after resampling\t", imgs_after_resamp.shape)
## Shape after resampling    (120, 140, 140)

3D Plotting

def make_mesh(image, threshold=-300, step_size=1):
    print ("Transposing surface")
    p = image.transpose(2,1,0)
    
    print ("Calculating surface")
    verts, faces, norm, val = measure.marching_cubes(p, threshold, step_size=step_size, allow_degenerate=True) 
    return verts, faces
def plotly_3d(verts, faces):
    x,y,z = zip(*verts) 
    
    print ("Drawing")
    
    # Make the colormap single color since the axes are positional not intensity. 
#    colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']
    colormap=['rgb(236, 236, 212)','rgb(236, 236, 212)']
    
    fig = plotly.figure_factory.create_trisurf(x=x,
                        y=y, 
                        z=z, 
                        plot_edges=False,
                        colormap=colormap,
                        simplices=faces,
                        backgroundcolor='rgb(64, 64, 64)',
                        title="Interactive Visualization")
    iplot(fig)
def plt_3d(verts, faces):
    print ("Drawing")
    x,y,z = zip(*verts) 
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)
    face_color = [1, 1, 0.9]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, max(x))
    ax.set_ylim(0, max(y))
    ax.set_zlim(0, max(z))
    ax.set_facecolor((0.7, 0.7, 0.7))
    plt.show()
    
v, f = make_mesh(imgs_after_resamp[::-1], 350)
## Transposing surface
## Calculating surface
plt_3d(v, f)
## Drawing

v, f = make_mesh(imgs_after_resamp[::-1], 350, 2)
## Transposing surface
## Calculating surface
plotly_3d(v, f)
## Drawing

Semi-Automatic Achilles Tendon Segmentation Using 3DSlicer

In order to train an automatic model based on a convolutional neural network we need a set of images already segmented. In this section there is a description of how to segmentate the Achilles tendon from an ankle CT using “a free, open source and multi-platform software package widely used for medical, biomedical, and related imaging research” called 3DSlicer. The following explanation is not aiming to serve as a manual, but to illustrate how semi-automatic segmentation can be performed. Therefore all the issues concerning the specific software way of using will not be specified.

First of all we need to load the data to 3DSlicer. In our case it is a set of 35 DICOM images pertaining to different axial DP slices of the ankle. Once they’ve been loaded, four different windows are displayed. The top-left one shows axial sections (horizontal slices, the images just as they are); the bottom-left one shows coronal sections (divides into ventral and dorsal); the bottom-right shows sagittal sections (divides into right and left), and the top-right builds a 3D volume of the ankle (enabling it at “Volume Rendering” module). At “Segment Editor” module we create three segments: one for the air (clay), one for the Achilles tendon (yellow) and another one for the rest of the ankle (green).

Threshold range coloured in clay

Threshold range coloured in clay

Using the “Paint” tool we plant some seeds for each segment on some slices for each section (4-5 images for each segment and section). Before that, we set a threshold range (from 2 to 1000, for example) to make the task easier since the out-of-range pixels won’t be considered. It’d be advisable to thoughtfully colour the “boundary slices” where the Achilles tendon joins the calcaneus bone.

Planting some seeds

Planting some seeds

After that, we use the “Grow from seeds” tool to classify every single pixel of every single slice in the threshold range into one of the three segments. We check the results and apply some modifications if necessary. We can view a provisional 3D construction of the tendon based on the segmentation so far by clicking on “Show 3D” and hiding the air and ankle segments.

Grown from seeds

Grown from seeds

As you can see, there are some pixels not belonging to the Achilles tendon coloured in yellow and some other belonging that are not coloured (little holes on the axial section yellow area). To correct this we remove the threshold range by disabling “Editable intensity range” and we use the “Erase” and “Smoothing” tools (besides “Paint”) on every axial section image. We choose the “Closing” smoothing method to fill the holes and the “Mean” smoothing method to smooth the edges. We can lower the segment’s colour opacity at the “Segmentations” module to help us hone the results.

Achilles tendon segmentation

Achilles tendon segmentation

We have finally segmented the Achilles tendon. Now we can remove the air and ankle segments and export the segmentation to a labelmap at the “Segmentations” module. 3D-Slicer also allows us to get the images’ radiomic features; to do it we switch to “Radiomics” module and generate the table using the labelmap.

Getting the radiomic features

Getting the radiomic features

Bibliography

[1]
P. Suetens, R. Verbeeck, D. Delaere, J. Nuyts, and B. Bijnens, “Model-based image segmentation: Methods and applications,” in AIME 91, 1991, pp. 3–24.
[2]
H. Chen, DICOM Processing and Segmentation in Python.” https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/
[3]
K. Ushey, J. Allaire, and Y. Tang, Reticulate: Interface to ’Python. 2022. Available: https://rstudio.github.io/reticulate/
[4]
A. C. Alamo, Como usar python en RStudio. 2020. Available: https://acca3003.medium.com/como-usar-python-en-rstudio-335f953c0efd
[5]
M. P. A. Starmans, S. R. van der Voort, J. M. Castillo Tovar, J. F. Veenland, S. Klein, and W. J. Niessen, “Chapter 18 - Radiomics: Data mining using quantitative medical image features,” in Handbook of medical image computing and computer assisted intervention, S. K. Zhou, D. Rueckert, and G. Fichtinger, Eds. Academic Press, 2020, pp. 429–456. doi: https://doi.org/10.1016/B978-0-12-816176-0.00023-5.
[6]
Farhana Sultana, Abu Sufian, and P. Dutta, “Evolution of Image Segmentation using Deep Convolutional Neural Network: A Survey,” 2020.
[7]
Shervin Minaee, Yuri Boykov, Fatih Porikli, Antonio Plaza, Nasser Kehtarnavaz, and D. Terzopoulos, Image Segmentation Using Deep Leanring: A Survey,” 2020.
[8]
K. K. D. Ramesh, G. Kiran Kumar, K. Swapna, Debabrata Datta, and S. S. Rajest, “A Review of Medical Image Segmentation Algorithms,” 2021.
[9]
Hyunseok Seo et al., “Machine Learning Techniques for Biomedical Image Segmentation: An Overview of Technical Aspects and Introduction to State-of-Art Applications (Author Manuscript),” 2020.